#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <sys/time.h>
#include <pthread.h>

/* 	Calculating the pi value using Monte Carlo simulation 	*/

long long n_iteractions;	// number of iteractions for each thread
long long * results;		// vector of results for each thread

/*	Funcion used by all threads 	*/
void * function ( void * args )
{
	// initialize variables
	int * args_d = (int *) args;
	int thread_id = *args_d;
	long long i, count = 0;

	// seed
	struct timeval time;
	struct drand48_data randBuffer;
	gettimeofday(&time, NULL);
	srand48_r(time.tv_usec * thread_id, &randBuffer);

	for(i = 0; i < n_iteractions; i++) {
		double x, y;

		drand48_r(&randBuffer, &x);
		drand48_r(&randBuffer, &y);
		x = 2 * x - 1;
		y = 2 * y - 1;

		//	condition to get one more point
		if(x * x + y * y <= 1)
			count++;
	}

	// returning the result of points inside the circle
	results[ thread_id ] = count;
	return NULL;
}

int main (int argc, char ** argv)
{	
	if(argc != 3) {
		printf("Usage: ./monte-carlo-par <number_of_iteractions> <number_of_threads> \n");
		exit(1);
	}

	/*	alloc and initialize variables	*/
	long long n = atoll(argv[1]);		// total number of iteractions
	int n_threads = atoll(argv[2]);		// number of threads
	n_iteractions = n / n_threads;		// number of iteractions for each thread

	// vector of threads and thread ids 
	pthread_t * threads = (pthread_t *) malloc( sizeof( pthread_t) * (n_threads + 1)); 
	int * thread_id = (int *) malloc( sizeof( int ) * (n_threads + 1) );

	// vector of results (for each thread save their result) 
	results = (long long *) malloc( sizeof( long long ) * (n_threads + 1) );		

	// counters
	long long count = 0;
	int  j;

	// create all threads
	for(j = 0; j < n_threads; j++) {
		thread_id [j] = j;
		if(pthread_create(&threads[j], NULL, function, &thread_id[j]))
			fprintf(stderr, "Error creating thread %d\n", j), exit(1);
	}

	// join all threads
	for(j = 0; j < n_threads; j++) {
		if(pthread_join(threads[j], NULL))
			fprintf(stderr, "Error joining thread %d\n", j), exit(1);

		count += results[j];
	}

	double pi = 4.0 * (double) count / (double) n ;
	printf("%.10lf\n", pi);

	// free memory
	free(threads);
	free(thread_id);
	free(results);

	return 0;
}
